XeF2 1 - 300eV Multi-orb processing¶
PH
v2b 27/10/20
- Updating plot styles to HV.
- Added branching ratios & more literature data.
- Added ion yield data.
v2 07/10/20 (updated and distributed 26/10/20)
- Running manually on Stimpy, ePSproc-v1.3-dev with multiJob class + RF updated experimental data analysis.
- (See N2O preliminary data analysis for template, http://localhost:8888/lab/tree/dev/ePSproc/classDev/ePSproc_multijob_class_tests_N2O_011020_Stimpy.ipynb)
v1 27/03/20
- Running manually on AntonJr, ePSprc-v1.2 env.
- Code adapted from recent multi-E chuncked job notebooks.
Setup¶
!hostname
!conda env list
import sys
import os
from pathlib import Path
import numpy as np
# import epsproc as ep
import xarray as xr
import matplotlib.pyplot as plt
from datetime import datetime as dt
timeString = dt.now()
# For module testing, include path to module here, otherwise use global installation
if sys.platform == "win32":
modPath = r'D:\code\github\ePSproc' # Win test machine
winFlag = True
else:
modPath = r'/home/femtolab/github/ePSproc/' # Linux test machine
winFlag = False
sys.path.append(modPath)
import epsproc as ep
# Plotters
from epsproc.plot import hvPlotters
# Multijob class dev code
from epsproc.classes.multiJob import ePSmultiJob
hvPlotters.setPlotters(width = 700)
# import bokeh
# import holoviews as hv
# hv.extension('bokeh')
Load data¶
# # Scan for subdirs, based on existing routine in getFiles()
fileBase = r'D:\projects\ePolyScat\xef2\XeF2_highRes_wf' # Data dir on Stimpy
data = ePSmultiJob(fileBase, verbose = 0)
keys = [4,5,6] # Set for 4d data only
data.scanFiles(keys = keys)
data.jobsSummary()
System properties¶
data.molSummary()
Orbitals 21 - 25 are the Xe(4d) core levels, with $A_{1g}$(SG), $E_{1g}$(PG) and $E_{2g}$(DG) components, where the classifcations are in full $D_{\infty h}$ and ePolyScat labels in parenthesis. The electronic structure was computed in Gamess (HF/SPK-QZP-rel), and converged to a bond-length of 1.9373A, and 4d energies as shown in the table above (~75 - 76eV). (See Molecular orbitals below for plots of the orbitals.)
Below: literature values. Comparing orb 21 to the lower 4d SO binding energy of 70.34 eV, the current computational results are off from the experimental binding energies by ~6.2 eV.
Source: Bancroft, G. M., P. A. Malmquist, S. Svensson, E. Basilier, U. Gelius, and K. Siegbahn. “Gas-Phase ESCA Studies of Valence and Core Levels in Xenon Difluoride and Xenon Tetrafluoride.” Inorganic Chemistry 17, no. 6 (June 1978): 1595–99. https://doi.org/10.1021/ic50184a040.
These values also compare well with the current experimental results - shown below from the file XeF2_4d_LV_115eV_binding_energy_extractions_with_atomic_Xe.png
Plot cross-sections and betas¶
These are from ePolyScat's getCro function, and are LF (unaligned ensemble) results. This provides a good, if general, overview. More detailed results - per orbital and symmetry - are given in the addendum. The structure/oscillations in some channels at higher energies remain to be explored/explained!
Overview (all symmetries, length gauge)¶
# Comparitive plot over datasets (all symmetries only)
Etype = 'Ehv' # Set for Eke or Ehv energy scale
Erange=[50, 300] # Plot range (full range if not passed to function below)
data.plotGetCroComp(pType='SIGMA', Etype = Etype, Erange = Erange)
# Comparative plot over datasets (all symmetries only)
data.plotGetCroComp(pType='BETA', Etype=Etype, Erange=Erange)
Branching ratios¶
# Stack XS data to new data structure
# NOTE: this is not quite correct, since it forces all data to same Ehv axis, but OK for a quick hack, and easy to pull out branching ratios
dsXS = xr.Dataset() # Set blank dataset, this is easier for stacking, probably
lText = []
for key in data.data.keys():
subset = data.data[key]['XS'].sel({'Sym':'All', 'XC':'SIGMA'}) # Set XS data, all syms only
# NOTE currently missing full dataset resolution for orb24, so try interp. (2eV not 1eV step size)
# Note dropna to ensure no NaNs, see http://xarray.pydata.org/en/latest/interpolation.html#interpolating-arrays-with-nan
if key == 'orb24':
subset = subset.dropna('Eke').interp(Eke = dsXS.Eke, method = 'cubic')
# dsXS[key] = subset.copy() # Set .copy() for safety here!
dsXS[data.data[key]['jobNotes']['orbLabel']] = subset.copy() # Set .copy() for safety here!
# Convert to Xarray & normalise
dsXS = dsXS.to_array().rename({'variable':'channel'}).squeeze()
dsXS['total'] = dsXS.sum('channel') # Sum over channels
# Normalise...
dsXS = dsXS/dsXS['total']
# Plot
Erange = [75, 250]
dsXS.swap_dims({'Eke':'Ehv'}).sel(**{Etype:slice(Erange[0], Erange[1])}, Type='L').plot.line(x='Ehv');
These are broadly similar over the different gauges. There's quite a bit of low-E structure, and things smooth out at higher-E, with ~2:1 ratio between SG:PG or DG channels.
TODO: detailed comparison with expt. Below are the results including SO states, which either need to be summed over (expt) or simulated (theory) for direct comparison! It looks like the $\Sigma:\Pi$ or $\Delta$ ratios are somewhat similar as is, maybe, some of the low-E structure.
Experimental data¶
- Set 1: RF processed results, 02/10/20
- Set 2: DH ion yields, 27/10/20
# Load experimental data
dataPathExpt = Path(r'D:\projects\XeF2_Soleil_2019\RF_data_analysis_021020')
# Set data attribs in dict similar to ePS results structure
dataFiles = {}
dataFiles['SIGMA'] = {'Data':'XS', 'File':r'XeF2_Xe4d_4d32_4d52_cross_sec_all_photon_energies_02102020.dat',
'States':['4d3/2', '4d5/2']}
dataFiles['BETA'] = {'Data':'Beta', 'File':r'XeF2_Xe4d_beta_all_photon_energies_02102020.dat',
'States':['\pi1/2 (4d3/2)', '\delta3/2 (4d3/2)', '\sigma1/2 (4d5/2)', '\pi3/2 (4d5/2)', '\delta5/2 (4d5/2)']}
# Update with ion data
dataFiles2 = {}
dataFiles2['ION-low'] = {'Data':'XS', 'File':r'XeF2_ion_yield_low_energy_cal.txt'}
dataFiles2['ION-high'] = {'Data':'XS', 'File':r'XeF2_ion_yield_high_energy_cal.txt'}
# Read data files and convert to Xarray
# 27/10/20 added quick hack to set 2nd array for ion yield data
dataList = []
dataList2 = []
for key in dataFiles:
dataIn = np.loadtxt(dataPathExpt/dataFiles[key]['File'])
# Convert to Xarray
dataXr = xr.DataArray(dataIn[:,1:], dims=('Ehv','State'), coords={'Ehv':dataIn[:,0], 'State':dataFiles[key]['States'][0:dataIn.shape[1]-1]})
dataXr.attrs = dataFiles[key]
dataXr.attrs['dataPath'] = dataPathExpt
dataXr.name = key
dataList.append(dataXr)
# Stack to Xarray
dataExpt = xr.concat(dataList, "XC")
dataExpt['XC'] = list(dataFiles.keys())
for key in dataFiles2:
dataIn = np.loadtxt(dataPathExpt/dataFiles2[key]['File'])
# Convert to Xarray
dataXr = xr.DataArray(dataIn[:,1].squeeze(), dims=('Ehv'), coords={'Ehv':dataIn[:,0]}) # Only 1D datasets in this case
dataXr.attrs = dataFiles2[key]
dataXr.attrs['dataPath'] = dataPathExpt
dataXr.name = key
dataList2.append(dataXr)
dataExpt2 = xr.concat(dataList2, "XC")
dataExpt2['XC'] = list(dataFiles2.keys())
# Quick plot
marker = 'x'
for item in dataExpt:
plt.figure() # Force new figure
item.dropna('State').plot.line(x='Ehv', marker=marker, ls=':') # Include dropna here to remove empty state dims
# Fix axis labels
if item.XC == 'SIGMA':
plt.ylabel('XS/Mb')
else:
plt.ylabel(r"$\beta_{LM}$")
Note dataset labels here - this dataset has XS summed over spin-orbit states (labelled 4d3/2 and 4d5/2), and betas for spin-orbit states.
# Quick plot - Ion yields
marker = 'x'
for item in dataExpt2:
plt.figure() # Force new figure
item.plot.line(x='Ehv') # Include dropna here to remove empty state dims
# Fix axis labels
plt.ylabel("Yield")
Comparison plots...¶
# Compare with computational results
pType = 'SIGMA'
Erange = [50, 250]
pltObj, lText = data.plotGetCroComp(pType=pType, Etype = Etype, Erange = Erange, returnHandles = True)
# Add expt data - cross-secions
scale = dataExpt.sel({'XC':pType}).max() / data.dataSets['orb21']['XS'].sel({'XC':pType, 'Type':'L'}).max() # Set scaling to match computational data
(dataExpt.sel({'XC':pType})/scale).dropna('State').plot.line(x='Ehv', marker=marker, ls=':');
# Add expt data - ion yields
ionData = 'ION-high'
scale = dataExpt2.sel({'XC':ionData}).max() / data.dataSets['orb21']['XS'].sel({'XC':pType, 'Type':'L'}).max() # Set scaling to match computational data
(dataExpt2.sel({'XC':ionData})/scale).plot.line(x='Ehv');
# Manual legend fix
lText.extend(dataExpt.sel({'XC':pType}).dropna('State').coords['State'].data)
lText.extend([ionData])
plt.legend(lText)
# Fix axis labels
if pType == 'SIGMA':
plt.ylabel('XS/Mb')
plt.title('XS (normalised to computational results)')
else:
plt.ylabel(r"$\beta_{LM}$")
plt.title(r"$\beta_{LM}$")
EDIT 27/10/20: now with ion yield data, which should be taken as the reliable experimental XS result. This doesn't show structure, but is broadly similar to the photoelectron measurements. Structure in photoelectron results to be ignored! Some of the previous comments still apply.
Original notes 13/10/20:
This is interesting... the structures in the computational XS look very similar to the experimental results, but on a much smaller energy scale. I suspect there could be many reasons for this in general (both computational and experimental)... but given that the betas match well (see below), this possibly suggests that this is a signature of multi-electron effects in the bound state.
Some reasons this might be the case:
- The ePS computations are single active electron type calculations, hence do not include any multi-electron/state-mixing effects.
- The good agreement with the betas suggests that the scattering potential, hence continuum wavefunction and angular overlap integrals, are, however, well modelled in the computations.
- This, to me, points to the radial overlap integrals in the computations being quite different to reality, hence anything which strongly influences the radial wavefunction of the bound and/or continuum state. (As a side-note, the bond length used in the calculations could also have an effect here, this was 1.937A, pretty close to the lit value of 1.988A; any other state-mixing, e.g. vibronic couplings, could also play a role of course, although these might generally be expected to affect the angular momentum composition too.)
- Given that the 4d states are mixed, and we also expect multi-electron effects, the bound state seems the most likely candidate for this.
Any suggestions/comments here would be appreciated, since this is rather speculative, and I don't have a good overview of the spectroscopy and general expectations here at all.
# Beta comparison plot over orbs
pType = 'BETA'
pltObj, lText = data.plotGetCroComp(pType=pType, Etype = Etype, Erange = Erange, returnHandles = True)
# Add expt data
dataExpt.sel({'XC':pType}).dropna('State').plot.line(x='Ehv', marker=marker, ls=':');
# Manual legend fix
lText.extend(dataExpt.sel({'XC':pType}).dropna('State').coords['State'].data)
plt.legend(lText)
# Fix axis labels
if pType == 'SIGMA':
plt.ylabel('XS/Mb')
plt.title('XS (normalised to computational results)')
else:
plt.ylabel(r"$\beta_{LM}$")
plt.title(r"$\beta_{LM}$")
As previously, applying an energy shift to the data gives a very nice agreement for the betas...
Eshift = 9.0 # NOTE - this breaks energy selection code later, for some reason - float type and/or tolerance in slice?
# Beta comparison plot over orbs
pType = 'BETA'
pltObj, lText = data.plotGetCroComp(pType=pType, Etype = Etype, Erange = Erange, Eshift = Eshift, returnHandles = True)
# Add expt data
dataExpt.sel({'XC':pType}).dropna('State').plot.line(x='Ehv', marker=marker, ls=':');
# Manual legend fix
lText.extend(dataExpt.sel({'XC':pType}).dropna('State').coords['State'].data)
plt.legend(lText)
# Fix axis labels
if pType == 'SIGMA':
plt.ylabel('XS/Mb')
else:
plt.ylabel(r"$\beta_{LM}$")
Addendum¶
Molecular orbitals¶
SEE: https://epsproc.readthedocs.io/en/dev/methods/ePSproc_orbPlot_tests_130520.html
from epsproc.vol import orbPlot
# chemPath = r'/home/femtolab/python/chem/tools/chemlab/chemlab' # Linux dev machine
chemPath = r'D:\temp\chemlab\chemlab' # Win dev machine
out = orbPlot.importChemLabQC(chemPath = chemPath)
# Test class
# filename = r"/home/femtolab/ePS/XeF2/electronic_structure/xef2_SPKrATZP_rel_geom.log" # XeF2 Gamess file OK
filename = r"D:\projects\ePolyScat\xef2\electronic_structure\xef2_SPKrATZP_rel_geom.log"
mo = orbPlot.molOrbPlotter(chemPath = chemPath, fileIn = filename)
# Set orbs to plot - currently handles list of strs only!
orbs = [21,22,23,24,25]
orbList = [str(item) for item in orbs]
[mo.calcOrb(orbN = orbN) for orbN in orbs];
# Plot orbs as independent figures
[mo.plotOrb(orbN = [item], interactive = False, showAtoms = False) for item in orbList];
This shows orbitals 21 - 25, the full set of 4d components. Note that (22,23) and (24,25) are degenerate pairs, and are treated as one orbital (with occ = 4) in ePolyScat, as noted previously.
Gauge results comparison¶
This shows a big change in absolute X-section between lenght (L) and velocity (V) gauge calculations (~order of magnitue), but the feature are appoximately the same in both cases (apologies for the poor plots here!).
data.plotGetCro(pType='SIGMA', Etype=Etype, Erange=Erange)
# Branching ratios
dsXS.swap_dims({'Eke':'Ehv'}).sel(**{Etype:slice(Erange[0], Erange[1])}).plot.line(x='Ehv', col='Type');
Beta parameters are pretty much idendical for both gauges, with a broad, shape-resonance type feature at lower photon energies, and lots of small-scale resonance structure to higher photon energies.
data.plotGetCro(pType='BETA', Etype=Etype, Erange=Erange)
Matrix elements¶
data.lmPlot(dataType = 'matE', Erange = Erange, Etype = Etype, thres = 0.1, logFlag = False, selDims = {'Type':'L'}, fillna = True)
Versions¶
import scooby
scooby.Report(additional=['epsproc', 'xarray', 'jupyter'])